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We introduce a dynamical model of coupled directed percolation systems with two particle species. The two 
species A and B are coupled asymmetrically in that A particles branch B particles whereas B particles prey on A 
particles. This model may describe epidemic spreading controlled by reactive immunization agents. We study 
nonequilibrium phase transitions with focused attention to the multicritical point where both species undergo 
the absorbing phase transition simultaneously. In one dimension, we find that the inhibitory coupling from B 
to A is irrelevant and the model belongs to the unidirectionally coupled directed percolation universality class. 
On the contrary, a mean field analysis predicts that the inhibitory coupling is relevant and a new universality 
appears with a variable dynamic exponent. Extensive numerical simulations on small-world networks confirm 
our predictions. 



PACS numbers: 64.60.Ht,05.70.Ln,05.40.-a,89.75.Da 

It is a challenging problem to establish the concept of the 
universality class in nonequilibrium phase transitions. In con- 
trast to equilibrium systems, nonequilibrium systems lack a 
general framework to begin with and display too diverse criti- 
cal phenomena to be classified into a few universality classes. 
Recently the absorbing phase transition has been attracting 
much theoretical interest since the universality class concept 
can be applied. It is a phase transition from an inactive (dead) 
phase into an active (live) phase, which may occur in systems 
with microscopic trapped (absorbing) states which the system, 
once trapped, cannot excape out of. 

Absorbing critical phenomena are categorized into a few 
universality classes depending on conservation in dynam- 
ics and symmetry between absorbing states (see for review 
Ref. li|]). The directed percolation (DP) class is the most 
well-established one, which encompasses a vast number of 
systems . The contact process |4] is an archetype of the 
DP class. In this model, each lattice site is either empty (0) 
or occupied by a particle (X). Each particle may annihilate 
spontaneously (X — > 0) or branch an offspring to a neighbor- 
ing site (X —> XX). The contact process was originally sug- 
gested as a model for epidemic spreading without immuniza- 
tion 1 4], where the particle represents an infected individual, 
who may recover spontaneously or infect its healthy neigh- 
bors. The vacuum state with all sites empty is the unique ab- 
sorbing state, which is one of the key ingredients in the DP 
class. 

Recently, coupled DP systems have been studied exten- 
sively In particular, the unidirectionally cou- 
pled DP (UCDP) process introduced by Tauber et al. |7|] was 
found to display a series of new multicritical phenomena. In 
the UCDP process, there are an infinite number of particle 
species Xk with k = 0, 1,2, • • • with dynamic rules character- 
ized by Xk — » 0, Xk — > XkXk, and X^ — > X^X^+i. The coupling 
is unidirectional (upward), linear, and excitatory. The species 
Xq is not affected by the others, hence its absorbing critical 



phenomena belong to the DP class. The other species X^q, 
however, are dressed by the critical fluctuations of particles of 
species Xy^ at the multicritical point where all species un- 
dergo the absorbing phase transition simultaneously. It turned 
out that their critical behaviors are described by the critical 
exponents distinct from those of the DP class 0, 0] . This 
novel critical behavior is destroyed and the DP scaling is re- 
covered as soon as we turn on the excitatory coupling in the 
other (downward) direction 001 . 

It is noteworthy that the UCDP process was proposed for 
describing an interesting roughening transition in an one- 
dimensional growth process |9], where the monomer-type 
restricted-solid-on-solid (RSOS) model is considered with the 
constraint that evaporation is allowed only at the edges (not 
terraces). It can be easily shown that the growth model is 
mapped to the UCDP process, except that the RSOS condition 
generates an extra inhibitory coupling in the downward direc- 
tion, i.e. the offspring production of X^ is suppressed by the 
presence of Xy for k' > k (see Sec.VA in |8]). The coupling 
should be at least bilinear, so the simplest process one may 
consider is X^X^+i — >Xk+\. Numerical studies |8, 9] revealed 
that the one-dimensional growth model displays the same crit- 
ical behavior as the UCDP, which may imply the irrelevance 
of the inhibitory coupling in coupled DP processes. 

In this Letter, we study a coupled DP system with two par- 
ticle species, which has a linear excitatory coupling in one 
direction and a bilinear inhibitory coupling in the other di- 
rection. We show that the inhibitory coupling is irrelevant 
only in low dimensions (UCDP class) but becomes relevant 
in high dimensions, at least where the mean field theory be- 
comes valid. Our model may also serve as a prototype model 
for epidemic spreading process in the presence of reactive im- 
munization agents, which will be explained shortly. 

We consider a coupled DP system with two particle species 
A and B on lattices or networks in general. Specifically we 
choose the contact process dynamics for both species in the 
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FIG. 1: A schematic phase diagram of the ACDP process in the 
(pa,Pb) plane with fixed A, and /j. 



noninteracting limit; each particle annihilates spontaneously 
or branches one offspring of the same kind to a neighboring 
site. The two species are coupled asymmetrically through 
branching and predation processes, that is, an A particle 
branches a B particle to a neighboring site (excitatory cou- 
pling), while a B particle preys on an A particle at the same 
site (inhibitory coupling). 

The dynamic rules are summarized as follows: (i) A (B) — ► 
with probability pa (pb), (ii) A — > AA with probability (1 — 
Pa)(1—X), (iii) B^BB with probability 1 - p B , (iv) A — > AB 
with probability (1 — pa)K and (v) AB — > B with probability p. 
The order parameters for the phase transition are the particle 
densities, Pa and p#. The model will be referred to as the 
asymmetrically coupled DP (ACDP) process. When p = 0, 
the coupling becomes unidirectional and the model reduces to 
the two-species case of the UCDP process. 

In the context of epidemiology, A particles correspond to 
antigens (or viruses) and B particles correspond to antibod- 
ies (or vaccines). Both antigens and antibodies can replicate 
themselves and may die out by themselves (processes (i), (ii), 
and (iii)). In addition, an antigen promotes an antibody (iv) 
and can be destroyed by the antibody (v). 

A schematic phase diagram of the ACDP process is shown 
in Fig.^ in the phase I where both annihilation probabilities 
are large, both species are inactive with pA = Pb = in the 
steady state. When ps is smaller, the B species becomes active 
with the A species remaining inactive. Hence, in the phase II, 
Pa = and pg ^ in the steady state. When pa is small, the A 
species is active, so is the B species regardless of the value of 
Pb due to the process A — > AB. So, both species are active with 
Pa ^ and p# in the phase III. Since Pa = throughout 
the phases I and II, the phase boundary line MR should be flat 
with ps = Pg which is a constant independent of pa, A,, and 
p, and is identical to the critical probability of the uncoupled 
system. 

The nature of the phase transitions is easily inferred. Across 
the line MR, the B species undergoes the absorbing transition 
with pA = 0. Across the line MP, the A species undergoes 
the transition while p# remains finite. Across the line MQ, 
both species undergoes the transition with p^ « Pa since the B 
particles are mostly slaved by the A particles through the A — > 



AB process for large ps region. Therefore, one may easily 
expect that the three critical lines should belong to the DP 
class. That was confirmed by numerical simulations, which 
we do not present here. Interesting critical phenomena can be 
observed along a line through the multicritical point M from 
the phase III to the phase I. In this Letter, we determine the 
location of M with high precision and mainly focus on density 
decay dynamics of both species at M, which may display a 
novel critical behavior distinct from the DP and also from the 
UCDP. 

First, we investigate the model using the mean field theory 
which assumes that local densities of A and B particles are 
uniform and uncorrected spatially with the mean values Pa 
and Pb, respectively. The mean field rate equations can be 
written as 

Pa = a A pA-pl-J>PAPB, 

Pb = flfiPfl-pi + ApA, (1) 

where the densities are rescaled to fix the coefficient of the 
p 2 terms to be unity and the parameters are functions of the 
model parameters. Higher order terms are neglected. 

The stationary state densities are obtained by solving 
Eq. with Pa = Pb — and hence the phase diagram is de- 
termined. It has the same structure as shown in Fig. \l\ The 
critical lines MR, MQ, and MP are given by a B — (a A < 0), 
«a = (as < 0), and ca = pas (as > 0), respectively. The 
multicritical point M locates at (oajOb) = (0,0). 

Approaching the multicritical point along a straight line 
(aA,aB) = (6,re) with fixed r inside the phase III, we find that 
the density vanishes as Pa,b ~ $ A B with the order parameter 
exponents Pa = 2 and p^ = 1. Spatial fluctuations in the order 
parameter are easily incorporated into the mean field theory 
by including diffusion terms V 2 pA,e in the mean field equa- 
tions. With those, one can show that the characteristic length 
scale diverges as ~ e~ V/1 B with the correlation length ex- 
ponent Va = Vb = 1 /2. The characteristic time scale diverges 
with the length scale as 1a,b ~ %a 'b w ^ tn tne dynamic expo- 
nent za=zb = 2. 

We note that the critical exponents for the B species are 
identical to the mean field DP exponents with p = 1,V = 1/2, 
and z = 2 (jl . On the other hand, the A species dynamics is 
dressed by the critical DP fluctuations and falls into the non- 
DP class. In the UCDP limit with p = 0, we have the opposite 
situation; A species falls into the DP class with Pa = 1 and 
B species into the non-DP class with Pb = 1/2 |7]. From the 
mean field analysis we conclude that the inhibitory coupling is 
relevant and that the ACDP forms an independent universality 
class distinct from the UCDP class. 

More interesting is the critical density decay dynamics at 
the multicritical point. Putting «a = og = in Eq. Q, after 
straightforward calculations, we obtain that the density decays 
algebraically in time as pA ~ t~ aA and p# ~ ?~ as with the 
exponent 0Cg = 1 and the variable exponent 

a A = max{2,p} . (2) 
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At /j — 2, logarithmic corrections appear as ~ 
t ~ 2 (2X In t ) ~ 1 and p B ~ t ~ 1 ( 1 + 1 / (2 In t ) ) . In the normal scal- 
ing, we expect that the density decay exponent is given by 
a = P/(zv) OL which does not hold in our model for /J > 2. 
This is another novel feature of the ACDP universality class. 

The mean field characteristics may be realized on regu- 
lar lattices in higher space dimensions than the upper critical 
dimension and also on various networks including complete 
graphs, random networks, and small-world networks I Kill . 
Recently, the small-world network is found to be useful and 
efficient in studying mean field critical phenomena 1 1 1 ] and it 
reflects some of network aspects in real biological and social 
systems. In this Letter, we adopt the small-world networks to 
test our mean-field result and also to get the information on 
the upper critical dimension of our model. 

The small- world network is characterized by the total num- 
ber of sites N, the average coordination number 2Z, and 
the so-called rewiring probability p r . It is constructed as 
follows fioll : First, N sites are arranged along an one- 
dimensional ring and each site is connected up to Z neighbor- 
ing sites. Second, each bond is rewired randomly with proba- 
bility p r . It interpolates the normal one-dimensional lattice at 
p r = and the random network at p r = 1 . In simulations, we 
choose Z = 10, p r = 0.5, and N <2x 10 6 . 

In networks, the concepts of dimensionality as well as 
length scale are meaningless. Botet et al. 11211 suggested that 
a correlation volume, instead of length, Z,v be a proper scaling 
variable and the criticality be described by the correlation vol- 
ume exponent V and the effective dynamic exponent z. More- 
over, they conjectured that the exponents are given by 



v = v 
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d u and z = z MF /d u 
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where v MF and z MF are the mean-field exponents and d u is 
the upper critical dimension. Although the conjecture has not 
been proved yet, it is shown to be valid for some equilibrium 
systems II UM2II . We will use this relation to find the upper 
critical dimension of our system. 

The phase boundary MR is given by a flat line of ps — p c B , 
which is the critical probability of the uncoupled system. 
We perform Monte Carlo simulations starting with the fully 
occupied configuration of B particles only. All numerical 
data shown below are obtained after average over, at least, 
several thousands of samples. As the particle density de- 
cays algebraically in time at criticality, we can easily locate 
p% = 0.48505(5) where p B - with a B = 0.98(5). This 
confirms that the critical line MR belongs to the mean-field 
DP class. 

Once p B being found, the multicritical point M can be lo- 
cated by examining the decay dynamics of p^, starting from 
the fully occupied configuration of A particles, as one varies 
Pa with ps = p% fixed. For the coupling constants, we fix 
A, = 0.5 and vary p . We first study the UCDP limit at p = 0. In 
Fig. 13(a), the particle density decay is presented at the multi- 
critical point at p c A = 0.3202(1 ). The density decay exponents 
are estimated as = 0.98(3) and Og = 0.49(2) as indicated 
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FIG. 2: (a) Density decay (N = 2 x 10 ) and (b) characteristic time 
scale at ^ = and p A = 0.3202. (c) Density decay (N = 2x 10 6 ) and 
(d) characteristic time scale at p = 1.0 and p A = 0.2902. 



by straight lines in the plot. Note that the results are consis- 
tent with the mean field results of the UCDP |7]. In Fig. 0(b), 
we present the scaling of the characteristic time x at which the 
survival probability becomes 1/2. It grows algebraically with 
the total number of sites N with the effective dynamic expo- 
nent za = 0.51(2) and zb = 0.52(2) as indicated by straight 
lines. Using Eq. © and the mean field dynamic exponents 
Za = zb = 2 [7], we obtain that the upper critical dimension- 
ality is given by d" = 4 for both species, which is consistent 
with the prediction for the UCDP class [7]. 

Turning on the inhibitory coupling (^i / 0), we find dif- 
ferent scaling exponents. In Fig.[2](c), the density decay is 
shown at the multicritical point p c A = 0.2902(2) at p = 1. As 
predicted in our mean field theory, the estimated decay expo- 
nent a B = 1.02(3) is consistent with the mean-field DP value 
and cia = 2.20(5) represents the non-DP class. We also esti- 
mate the effective dynamic exponents za = 0.33(1) ~ 1/3 and 
z B = 0.50(1) in Fig. 0(d). Since z MF = 2 for both species, 
Eq. (|3} suggests that the two species have different upper crit- 
ical dimensions, i.e. d\ = 6 and d\ = 4. It is remarkable to 
see that the two interacting species begin to show the mean 
field critical behaviors at different dimensionality. A sim- 
ple power counting for coupled Langevin equations obtained 
from Eq. ([0 by adding diffusion terms and DP-type multi- 
plicative noises reproduces the above results, but more efforts 
are needed to understand this peculiar property both analyti- 
cally and numerically. 

We performed the similar analysis for various values of p. 
The results for the decay exponents and the effective dynamic 
exponents are presented in Fig. [3] The numerical results are in 
good agreement with the mean field results. Near/j = 0, e.g., 
at p = 0.25, the exponents seem to deviate from the mean field 
values. These are attributed to crossovers due to the nearby 
UCDP fixed point. The mean field theory and the numerical 
results suggest that the ACDP process with p ^ form a dis- 
tinct universality class different from the UCDP class. 

The previous study poses a question on the relevance of the 
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FIG. 3: (a) The density decay exponent and (b) the effective dynamic 
exponents at the multicritical points. The solid (dashed) lines are 
guides to eyes and indicate the mean field results for the A (B) species 
of the ACDP process. The solid (dashed) arrows point toward the 
corresponding exponent values of the mean field UCDP class. 




FIG. 4: (a) Density decay at the multicritical point of the ACDP 
model with X = 0.2 and /j = 0.2 on the one-dimensional lattice with 
N = 10 sites, (b) Characteristic time scales at the multicritical point. 



inhibitory coupling at lower dimensions H. We check the 
multicritical behavior of our model on the one-dimensional 
lattice. We fix the coupling constants X = 0.2 and /j = 0.2. 
Numerical simulations show that the multicritical point is lo- 
cated at (j>a,Pb) = (0.14627 (2), 0.23267). The density de- 
cay and the characteristic time scales are presented in Fig. [3] 
The decay exponents are estimated as oca = 0.159(5) and 
Ob = 0.10(2), which are consistent with the results for the 
UCDP class in one dimension |7]. The dynamic exponents 
are esimated as za—Zb — 1-59(2), which are consistent with 
the DP value. These results suggest that the ACDP process 
belongs to the UCDP class in one dimension. That is, the in- 
hibitory coupling is irrelevant in one dimension in the contrary 
to the higher dimensional cases. 



In summary, we have introduced the ACDP process of two 
particle species and investigated its multicritical behaviors. In 
the mean field theory, we show that the inhibitory coupling 
is relevant and the ACDP process forms the different univer- 
sality class distinct from the UCDP process: The B species 
displays the DP type critical behavior while the A species the 
dressed non-DP type critical behavior. Furthermore, the criti- 
cal dynamics of the A species is anomalous with the variable 
critical exponent. The mean field results are confirmed by the 
numerical simulation studies of the ACDP on the small-world 
networks. Using the conjecture in Eq. Q, we also find that 
the upper critical dimensionality for each species is given by 
d\ = 6 and dg = 4. On the other hand, the numerical simu- 
lation studies in one dimension show that the inhibitory cou- 
pling is irrelevant and the ACDP process belongs to the UCDP 
class: The A species displays the DP type critical behavior 
while the B species the dressed non-DP type critical behavior. 

For future works, we suggest that it would be interesting 
to study the model in higher dimensional lattices to find the 
dimensionality at which the inhibitory coupling becomes rel- 
evant and to cofirm the upper critical dimensionalities. The 
multi-species generalization would be also interesting fl3ll . 
As an application, the phase transition of the ACDP process 
on general complex networks may yield richer critical dynam- 
ics. 
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